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We have investigated the effect on phonon energy transport in mesoscopic systems and the reduc- 
tion in the thermal conductance in the quantum limit due to phonon scattering by surface roughness 
using full 3-dimensional elasticity theory for an elastic beam with a rectangular cross-section. At 
low frequencies we find power laws for the scattering coefficients that are strongly mode dependent, 
^ and different from the dependence, deriving from Rayleigh scattering of scalar waves, that is 

'^'h , often assumed. The scattering gives contributions to the reduction in thermal conductance with 

the same power laws. At higher frequencies the scattering coefficients becomes large at the onset 
' frequency of each mode due to the flat dispersion here. We use our results to attempt a quantitative 

understanding of the suppression of the thermal conductance from the universal value observed in 
experiment. 

' PACS numbers; 63.22.-|-m, 63.50.-|-x, 68.65.-k, 43.20.Fn 

(p : I. INTRODUCTION 

Landauer's formulation of quantum transport showed that when elastic spattering dominates, the electrical con- 
^ ' ductance can be related to the transmission coefficient of the electron waveail. In the ideal case of no scattering, this 
^ [ leads to a universal conductance that is quantized in units oi e^/h at low temperatures, with an additional quantum 
of conductance added as each channel or mode of the conductance pathway opens up. The application-af_aimilar ideas 
to the phonon counterpart, namely thermal conductance, was recently derived by a number of authorsoLlu, and is now 
recognizedu to be related to earlier work on the entropy transport at low temperaturesQ. Rego and KirczenowJriave 
^ ' extended the concept of the universality of the thermal conductance to particles of arbitrary statistics (anyons)Ll. 

In the case of electrical resistance, the chemical potential or the number of conducting modes can be varied at 
I ■ very low temperatures, giving sharp jumps between various quantized values of the resistance. On the other hand, 
^ , thermal transport by phonons necessarily requires nonzero temperatures to populate the modes of the conducting 

■ pathway, and the width of the Bosc distribution function smears out the quantization of the conductance. Only at 
very low temperatures, where just the modes of the conducting pathway with zero frequency at long wavelengths 
contribute to the thermal conductance, the quantization of the ideal conductance becomes apparent in a universal 

■ thermal conductance NqKu with = {tt^ / ?>)k^T / h the universal conductance per mode, with fcs Boltzmann's 
. constant and h Planck's constant, and A'o is the number of modes with zero frequency at long wavelengths, which 
' is four for a freely suspended elastic beam connecting the two thermal reservoirs. Note that this value of the low 

■~>^ temperature conductance in the absence of scattering is independent of the dimensions and elastic properties of the 
thermal pathway. 

Ch a low temperature thermal conductance consistent with the predicted universal value was measured by Schwab et 

I al.El in experiments on a lithographically defined mesoscopic suspended beam (of dimensions about 1/im x 200nm x 
■ O 60nni). Whilst their elegant experiment displays the universality of ballistic phonon transport, the experiment also 
C showed a decrease in the thermal conductance below the universal value in the temperature range of 0.08K < T < 0.4K 
that cannot be explained by the ballistic theory, since in this theory an increase in the thermal conductance is expected 
as the temperature is raised and more modes are excited. The decrease in thermal conductance is presumably 
associated with the scattering of the thermal phonons, and can be understood using the ideas of Landauer in terms 
of the scattering coefficient of the vibrational waves. This is the topic of the present paper. 
1 In this paper we calculate the effect on the low temperature thermal conductance of the scattering of the thermal 
' phonons by surface roughness, which is likely to be the major source of scattering in mesoscopic samples. The scatter- 
ing of scalar waves, described by the simple wave equation, in waveguides with routda, siirface s has been investigated 
by many workers, including ourselves, using both numerical and analytic methodsfli3't2ltiHl3. However, for the low 
frequency modes of interest in the low temperature thermal conductance, the physical vibrational waves have quite 
different properties than the waves in the scalar model. For example the dispersion relations of the modes are dif- 
ferent, with two of the four modes with zero long- wavelength frequency having a quadratic dispersion at small wave 
vectors, rather than the linear dependence given by the simple scalar theory. To understand the experimental results 
quantitatively, a more accurate treatment of the vibrational waves is needed. At low temperatures the wavelengths of 
the thermally excited modes are large compared with the atomic spacing, and so a treatment based on the equations 
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of macroscopic elasticity theory is appropriate. Blencowelli has considered the scattering of elastic waves in a thin 
plate waveguide with rough surfaces, but prior to our work, the scattering of elastic waves confined in a beam-like 
wave guide with rough surfaces has not been considered. 

Previously, we havCjiiivestigated the effect of surface scattering on the low temperature thermal conductance using 
the scalar wave modelEHI. In that paper we noted the apparent discrepancy between the results of scalar model with 
a simple assumption for the nature of the surface roughness and the data by Schwab et al. below a temperature of 
O.IK: the data seemed to show a delay of the onset of scattering as the temperature increased that was not predicted 
by the model. However, since the scalar model does not properly account for the properties of the elastic waves, it 
was not clear whether this discrepancy is due to an inadequate modelling of the surface roughness, or the flaw in the 
description of the waves themselves. To resolve this matter, and obtain a more accurate account of the scattering of 
the waves by rough surfaces, we develop a theory of based on the full elasticity equations, and use this tp-icalculate 
the thermal conductance at low temperatures. A short version of this work has been previously publishedtJ. 

In Section II, the scattering of elastic waves confined to a beam of rectangular cross-section with rough surfaces 
is calculated using the full three dimensional elasticity theory. We use a Green theorem approach, and calculate the 
scattering coefficient to quadratic order in the amplitude of the surface roughness. These results are quite general, 
but rather intractable for further progress, since the structure of the modes in an elastic beam cannot be determined 
in closed form. Thus in Section [II we reduce the expressions to a thin plate limit to provide a closed form for 
the displacement fields, and to obtain analytical expressions for the scattering behavior. In Section IV the general 
behavior of the scattering and the effect on the thermal conductance is analyzed in detail, using a simple description of 
the surface roughness, to investigate the physical consequences of the novel features of the elastic waves. In Section ^ 
we use our theory to attempt to fit the data of Schwab et alB using more realistic descriptions of the surface roughness. 
A number of the more difficult issues that arise in the elasticity theory are described in appendices. 

Although our main interest is the scattering of thermally excited vibrational waves in mesoscopic systems at low 
temperatures, the formulation of the surface scattering is quite general, and can be applied to other situations, such 
as the scattering of mechanically excited modes in macroscopic samples for example. 



II. GENERAL FORMALISM 



A. The model 



The main focus of this paper is the effect of surface roughness on the low temperature thermal conductance of 
mesoscopic structures. The geometry we consider is a freely suspended elastic beam, which we call the bridge, 
connecting two thermal reservoirs. We will consider a beam of rectangular cross-section of dimensions width W (in 
the y direction) and depth d (in the z direction). Mesoscopic structures are often produced lithographically from 
epitaxially grown material. We choose a convention that the depth is the dimension in the growth direction, and the 
width in the lithographically defined transverse direction. We define the length of the rectangular beam of nominally 
uniform cross section as L. In practice the bridge may be joined to the reservoirs smoothly, by a portion of continuously 
growing width, to eliminate or reduce the scattering of the vibration modes off a sharp junction. We will suppose 
that the scattering by roughness is important only in ^Qjpe narrower portion of length L. 

The thermal conductance is given by the expressionB'B'Q 



where ujm is the cutoff frequency of the m th mode, /3 — l/{kBT), and T is the temperature. The integration is over 
the frequency ui of the modes m that propagate in the structure. The transmission coefficient is unity for the ideal 
case. Any scattering reduces the thermal conductance, and scattering of the lowest modes can reduce the conductance 
below the universal value at low temperatures. 

To actually perform the scattering calculation we imbed the rough beam of length L in an infinite beam of the same 
cross section but with smooth surfaces outside of the region of length L, Fig. |l|. Thus the mathematical calculation is 
the scattering of a wave incident from x — — oo on a rough portion of the beam with surfaces a.t y = ±W/2 ± /i {x, z) 
and at z — ±c?/2 ± /2 {x,y), with the height functions /i,2, defining the roughness, nonzero only in a finite region 
< a; < L. Forward scattering is evaluated from the intensity of waves as a; ^ -l-oo, and backward scattering from 
the intensity of waves as a; — > — oo. 

To calculate the scattering amplitude, we take a Green function approach similar to our previous work on the scalar 
mode]ll3. 
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FIG. 1: Top: Three dimensional elastic beam with rectangular cross-section. The rough surfaces are on the top, bottom, 
sides. Bottom: Side view of the mathematical model of the structure actually used for the scattering calculation. 
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The displacement field u away from any sources satisfies the wave equation: 

pd^^u, = a,T„ (2) 

where p is the mass density, and 

Tij = CijkidkUi (3) 

is the stress tensor field with Cijki the elastic modulus tensor. The subscript i runs over the 3 Cartesian coordinates, we 
use the symbol to denote the derivative d/dx etc., and repeated indices are to be summed over. The displacement 
field satisfies stress free boundary conditions at the surfaces 

T,,n,\s = (4) 

where S denotes the surface boundaries and nj is normal to the surface. Assuming harmonic time dependence at 
frequency uj, Eq. becomes 

puj'^ui + Cijkidjdkui = 0. (5) 

We approximate the material of the system as an isotropic solid. Then the elastic modulus tensor is 

Cijki = ^StjSki + p (StkSji + 6ii6kj) (6) 

where A and p are Lame constants (p is also the shear modulus) 

Ea/ {I + <7){l -2(7), p=^E/2{l + (j) (7) 

with E Young's modulus and a the Poisson ratio. 

Even in a rectangular beam geometry the displacement fields in the propagating modes yielded by these equations 
are complicated, and cannot be found analytically. The modes can be grouped into four classes according to their 
signature under the parity operations y —y and z —z. Some modes show regions of anomalous dispersion where 
the group velocity duj/dk is negative: these regions require a careful examination of the notions of "forward" and 
"backward" scattering for the waves. The lowest frequency mode of each class has a frequency that tends to zero 
at small wave number. These four modes are the only ones excited at low enough temperature, and are the ones 
contributing to the universal thermal conductance. The structure of these modes at small wave numbers is simple 
and can be calculated using familiar macroscopic arguments of elasticity theory: they are the compression, torsion, 
and (two orthogonal) bending modes. 
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We define a Green function Giq{'x, x'; t, t') to satisfy the wave equation with a source term —5iq5{yi — x')S{t — t'), 
and ^ijq to be the corresponding stress 

^ijq = CijkldkGlq. (8) 

It is convenient to introduce the frequency space version of the Green function 

G.,(x;x';i,0- / ^G.,(x; x'; c^)e-'-(*-*') , (9) 

with a similar expression defining ryq(x, x'; w). Inserting G, F, and the source term into Eq. gives 

ptj^G,;,(x, x'; tj) + (9jTy,(x, x'; tj) = -6iq6 (x - x') (10) 

where x is the observation coordinate and x' is the source coordinate. 

Equations (^ and (|l0| ) lead to Green's theorem expressing the displacement field at frequency u in terms of a 
surface integral 

ii,(x) = / K-T,, (x') G,q (x', x; t^) - n'.u, (x') T,,, (x', x; to)] dS' . (11) 

We are free to choose any closed integration surface S' . One choice is to use the physical rough surface thereby 
eliminating the first term in Eq. ( [Ill ) due to the boundary condition Eq. (jj). However, the resulting integration over 
the rough surface is not easy. Instead, we integrate over the smoothed surfaces at y = ztW/2 and z = zLd/2 and 
impose the boundary conditions on the Green function to be stress free on these smoothed surfaces 

r,,,n,|^ = 0, (12) 

together with cross sections at x' ±oo to close the surface. 

The total field u can be written as the sum of incident and scattered waves 

u = u'" + u'". (13) 

It can be shown (see Appendix ^) that the integration over the sections at x' ztoo on the right hand side of Eq. 
( pi] ) just gives u^". In the integration over the smoothed surfaces at y = ±W/2 and z = ±d/2 the second term in the 
integrand vanishes due to Eq. (|2|). Thus we find the expression for the scattered field 



<(x) = / [n;ry (x') G,q (x', x; ^)] dS' , (14) 
Js 

with the surface S the smoothed surfaces y = ztW/2 and z = ztd/2. The stress field on the smoothed surface is 
evaluated by expanding about its value on the rough surfaces, where Eq. (|^) applies. 

The rest of the section goes as follows: firstly, we find an explicit expression for the Green function with stress free 
boundary conditions; then we apply the boundary perturbation method to project the stress at the rough surfaces 
onto the smooth surfaces by expanding the stress free boundary terms around the smooth surfaces using the small 
roughness as the expansion parameter; and finally we evaluate the strength of the scattered waves to give the scattering 
coefficient. 



B. Green function 

We evaluate Gig(x, x'; w) as an expansion in the complete orthonormal set of normal modes u'^'^'™^ (x) in the ideal 
geometry, which satisfy Eq. (|^) and stress free boundaries at the smooth surfaces. Here k is the wave number in the 
X direction, and m labels the branch of the dispersion curve. We define ujmik) as the frequency of the mode m at 
wave number k in the ideal geometry. The modes satisfy the completeness relation 

J2 I ^uf-)(x')*4'^^'")(x) =5,,<5(x-x'). (15) 
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Substituting this expression on the right hand side of Eq. (|10[) leads to the expression for the Green function 



Gi,(x',x;cj) = 2^ 



dk 



{k.7n) 



{uj + ie) -ujl^ik) 



(16) 



where we write 



,ikx 



(17) 



with 



,(fe,m) 



giving the transverse dependence of the displacement field. In Eq. ( [T6| ) e is a positive infinitesimal number 
to incorporate causality, Giq{x,x' ;t,t') = for t < t' . 

Equation ( p^ ) can now be evaluated by contour integration. The integrand has poles labelled by an index n near 
values fc = fc„ on the real axis which are given by solutions to the dispersion relation tOmikn) = lu for all branches to. 
(We take an incident wave with lu > 0.) Note that for branches with regions of anomalous dispersion there may be 
more than one solution to this equation for some to, so that the index n is not identical to the branch index m. The 
poles are shifted slightly off the real axis by the infinitesimal e in Eq. (119), and are given by expanding about kn 



with Ug"-* the group velocity at the nth pole dujm/dk\i^^i^ . Notice the poles are in the upper half plane for Vg'^' > 0, 

(n) 

and in the lower half plane for < 0. 

Now we can perform the k integration by complex integration. Consider first the case, x > x' . The contour must 
be closed in the upper half plane so that the contribution from the semicircle at large |fc| vanishes. The contour 

integration then picks up contributions from the poles in the upper half plane, i.e. wave numbers with v'g^ > 0. On 
the other hand, for x < x' the contour must be closed in the lower half plane and it is poles at wave numbers with 

in) 

Wg < that give nonzero residue. Forward scattering or backscattering is thus seen to be determined by the sign of 

fn) 

the group velocity Vg rather than by the sign of as indeed would be expected physically. 
Evaluating the residues gives the expression for the Green function: 



(«) 



G,,g(x',x;w) = i 



^i")(xO*4"^(x) 

o (") ' 



(18) 



where u^f'\x) is written for u^'"^\:x.) at the value of the wave number k ~ kn satisfying ajmikn) = ^- The prime on 



the sum is used to denote the fact that the sum runs over intersections 



th vf' > for X > 



and over n with 



Va^^ < for X < x'. 



(n) 

The group velocity Vg does not have an analytical expression for a rectangular beam, and is obtained by differ- 
entiating the dispersion curve which must be found numerically. Alternatively, to avoid numerical differentiation, we 



can rewrite Wg"' in terms of the average power flow in mode n. Since u\'^' is normalized, the power P„ in mode n can 
be written as 



,(«) 



■Re 



-ia;Ti"^M-"^* I dydz 



2„(«) 



(19) 



the first expression of the equality expressing the energy flux in terms of the rate of work done across a section, and 
the second in terms of the group velocity and the average energy density evaluated as twice the average kinetic energy. 

(n) 

Then Vg can be evaluated in terms of P„ as 



and Pn has an expression directly in terms of displacement field given by the first equality in Eq. ( [l9| ) 

1 



Re 



dydz 



(20) 



(21) 



This expression for Ug"^ can also be derived directly from the equations of motiorJi; 
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C. Boundary perturbation 

In this section we show the boundary perturbation technique for the rough surfaces on the sides (i.e. the xz 
boundary planes). We work out the scattering coefficient expUcitly for the surface near y = W/2. The surface near 
y = —W/2 wiU give a similar contribution and, assuming uncorrelated roughness on the two surfaces, is accounted 
for by multiplying the single-surface scattering rate by 2 at the end of the calculation. The results for the top and 
bottom surfaces can be obtained by interchanging y and z whenever they occur in the indices in the displacement 
fields and stress tensors in the calculation below. 

In order to calculate the stress on the smooth surface appearing in Eq. (|l^, we expand the stress in a Taylor 
series about the flat surface, and impose stress free boundary conditions at the rough surface which is the small 
distance /i away. We also assume fi is differentiable. 

The unit vector n normal to the rough boundaries to first order in /i is 

n~y-dxfi{x,z)x~dzfi{x,z)z. (22) 
Then the stress free surface boundary conditions Eq. (^ can be written 

[Tiy - dxfi [x, z) - d^fi (x, z) Tiz]y^wj^f = (23) 
Now expanding Eq. ( p3|) in the neighborhood of ?; = W/2 and taking only the lowest order in /i and /{, we obtain 

Txyly^w — [dxfi {x, z) + dzfi (x, z) - fi {x, z) dyTr^y)\y^w (24) 
Tzyly^w ~ idJi{x,z)Tz, + dJiix,z)Tzz- fi{x,z)dyTzy)\y^„ (25) 

Tyy\y=}v ^ -fl (x, z) dyTyyly^w (26) 

where the first two expressions for Txy and Tzy have been used to simplify Tyy. Since the terms on the right hand 
side of Eqs. (p^p6|) are explicitly first order in the small parameter /i, the stress field on the right hand side can 
be evaluated at zeroth order, i.e. for ideal smooth surfaces. These results are used in Eq. (H). 



D. Scattering coefficient 



We now evaluate the expression for the scattered field given by an integration over the beam surfaces Eq. (p^. To 
calculate the scattering coefficient, we consider an incident wave of unit amplitude in a single mode m. Again in this 
section we will outline the calculation for the scattering by the single surface at y = W/2, and will include the effects 
of the other surfaces at the end. We therefore have 



K (x) 



[Tty (x') Giq (x', x; uj)] dx'dz'. 



(27) 



We can now evaluate the forward and backscattering amplitudes by using Eq. (18) for the Green function in Eq. 
(p7|), and evaluating the scattered wave at large positive and negative x 



{x oo, y, z) 



dx' 



dz 



2pUJ Vg"^ 



T,y{^)uT^ (X') 



,(») 



V - — 



(x) 



(28) 



(x^ -oo,?;, z) 



dx' 



dz' 



E 



r, (n) 

ZpU) Vg 



T,y{^')u'r^ (X') 



(x). 



(29) 



The stress tensor corresponding to the full displacement field of the wave is evaluated from Eqs. ( p4[]2q ). Since 
these expressions explicitly include the small roughness amplitude /i on the right hand side, to calculate the scattering 
at lowest order in the roughness amplitude it is sufficient to replace all Ty on the right hand side by the value T'/^™'' 
in the incident mode m. From Eqs. (p8|), ( p9| ) we see that u'''^(x) is expressed as a sum over modes u(")(x), and the 
coefficient of each mode is then the scattering amplitude t„^m from incident mode m into mode n, so that 



dz 



2pu; -yg"-* 



(x)ur (x) 



(n) 



y=W/2 



(30) 
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where we can now extend the integration Hmit to ±00 since /i, and so the integrand, is zero outside the domain of 
roughness < a; < L. Again mode indices n for which Vg""^ > represent the forward-scattered waves and those with 

(n) 

Vg < the backward-scattered waves. 

Now use the expression for the stress tensor on the smooth surfaces obtained in the previous section Eqs. ( p^^2^ ) 
and integrate the resulting expressions by parts with respect to a; or z to rewrite the terms in dxfi and dzfi as 
integrations over /i. After these manipulations we find tn^m can be written 



dx / dzfi (a;,z)r(™'") {x,z) 



(31) 



where 



Applying the equations of motion Eq. (||) and remembering 
somewhat simpler expression 



(m) 



y=w/2 



(32) 

V=W/2 

for all i and for all x,z leads to the 



r(™'")(a;,z) = [-pcj2 (4")^")* 



(33) 



Notice that the scattering separates into a kinetic term (the first line) and a stress term (the second line). 

The above form for r*-™'") is still neither instructive nor practical for numerical evaluation. It can be further 
simplified using the expressions Eqs. (^ and (Q) for the stress tensor in terms of the displacements. First we use the 
boundary condition Tyy'^ = for the y-stress to give aX y — W/2 

This can be used to simplify the expressions for the other components of the stress tensors at y = W/2 

E 



rp(m) 

XX 



rp{m) ^ 



-i(m) 



(1 


-cr^) 






(1 






E 



2(1 +(7) 



(9,4") + (79,4™)), 
(a9,4™)+9,4™)), 

(9,4") + 9.4™)). 



Inverting these gives at y = W/2 



E 



9.4"^ = i(Tir)-cTri;")), 



E 

a + a 7/(™) - ^iiizlr^™) 

'-'x^z ^ i-'z"-x ~ ^xz ■ 



(34) 

(35) 
(36) 
(37) 

(38) 
(39) 
(40) 



We emphasize that Eqs. (| 
the material. 

Using these results we get 



are only true on the stress free boundaries, and are not generally true in the bulk of 



2puj Wg") 



dz/i(fc„-A:„,z)f(™'")(z) 



(41) 
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with 



1 



1 



(rf{m)rp(n)* , rp{m)rp(n)*\ _ ( rp{m,)rp(n)* I rTi(m)/Ti(n)*\ yi(m)^(ri)* 

V a;a; js2: ^ zz zz } \ zz xx ' a^a; zz / xz zx 

(")/ 



where we have introduced the exphcit x dependence of u]"^ (x) as in Eq. and the stress tensor 

7^/"^(x)=T.,(y,z)e''=-, 



(42) 



(43) 



so that the x' integration is just the Fourier transform / of the roughness function, and F is a function of the z 
coordinate only. 

Alternatively, using Eqs. (|3^^7|) we can derive an expression explicitly in the displacement fields, which is useful 
for numerical evaluation, 



p(m,n) _ I 



2 / ifm) i(n)* I i(m) i(n) 

Wx 4>x +4>l '4>y ' 



2fi 



(l-a) 



(44) 



Tha^cattering rate is given by multiplying \tn,m\ by the ratio of the group velocities in the scattered and incident 
wavesEj. We also treat the roughness of the surface statistically, and take an ensemble average (denoted by angular 
brackets) to give the final expression for the scattering rate "fn,m from mode m to mode n by the per unit length of 
single rough surface at y = W/2 given by 



(n) 

|2\ 



1 



4p2 



9 (m) (n) 



(45) 



We are interested in the reduction of the phonon heat transport due to rough surfaces. Only the backscattered 
waves (those with Ug"^ < 0) reduce the amount of heat transmitted. Thus we define 7m,, the thermal attenuation 
coefficient of mode m per unit length, to be the sum of the scattering rates from the incident mode m to all possible 
backscattered modes, per unit length of rough surface. This can be written for scattering off the single rough surface 
considered so far 



E 

n 

t>i,"'<0 



1 



A 2 2 ("^ 



dz/i(fc„~A;„,z)f("-")(z) 



(46) 



To include the second rough side surface, assuming uncorrelated roughness, we simply have to multiply the expres- 
sion for 7m by a factor of 2. The expression for scattering off the top and bottom surfaces, if these are rough too, 
can be derived in a similar manner, and the result may be obtained by exchanging y and z in Eq. (^). The total 
scattering rate is the sum of the scattering off all the surfaces. 

We have assumed that the amplitude of the surface roughness is small, allowing us to use perturbation theory to 
derive the above expressions. In this case, the attenuation coefficient gives an exponential decay of the energy flux, 
so that the energy flux transmission coefficient is 



Tm = exp [-7mi] 



(47) 



III. THIN PLATE LIMIT 



Although the expression in the previous section is general and applicable to any rectangular waveguide with rough 
surfaces, there are no closed-form expressions for the displacement fields in general, and so a direct evaluation of 
the scattering would have to be done completely numerically. Here, we instead use the thin plate approximation 
d <C VFt3E3, which yields closed form expressions for the displacement fields of the modes (in terms of a dispersion 
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curves Umik) given by numerical solution of a simple transcendental equation). The thin plate limit captures the 
important properties of the elastic modes, for example the quadratic dispersion of the bending modes at small wave 
numbers, and regions of anomalous dispersion, as well as providing analytical expressions enabling us to do further 
analysis of the scattering. The thin plate theory is applicable where the thickness of the sample is much less than the 
width and the wavelengths are much greater than the thickness, which is the case for many mesoscopic systems at 
low temperatures. l— , 

The use of the thin plate limit for mesoscopic structures was proposed in referencetS, where the calculation of the 
structure of the modes is described in more detail. It is found that the modes can be separated into two classes: 
in-plane modes, where the polarization of the displacement is largely in the xy plane (together with small strains in 
the z-direction given by the Poisson effect) and the displacement field is completely specified by giving the vertically 
averaged horizontal displacement components Ux{x,y) and Uy{x,y); and flexural modes, where the displacement is 
primarily in the z direction and is specified by a vertical displacement field Uz{x,y). Within each class we can 
further distinguish the modes by their parity under y — > —y. For the in-plane modes we define the mode as even if 
Ux {x, —y) — Ux [x, y) and odd if Ux {x, —y) — —Ux {x, y). Similarly, the even flexural modes have Uz [x, —y) — Uz {x, y) 
and the odd modes have Uz {x,—y) = —Uz {x,y). As in the general case, there are four branches of the dispersion 
curves that tend to zero frequency as the wave number goes to zero, corresponding to one mode from each of these 
classes. The low frequency, even in plane mode corresponds to the compression mode, and the odd mode to a bending 
mode. The low frequency even, flexural mode corresponds to the second bending mode, and the low frequency odd 
flexural mode is the torsion mode. l— . 

Explicit expressions for the displacement fields can be obtained using the method described in referencellS. For the 
in-plane modes we find, up to a normalization factor Ai that is common to both even and odd parity waves, the even 
modes 



Ux {x,y) = ikAi 



k^-xl (X2W\ fxiW 



Akx 



(48) 



Uy {x,y) = Ai 



2 ^,2 



— cos — — sm (xi 



and the odd modes 



Ux {x,y) = ikAi 



2x1 V 2 ; 

k'-xl 



iXiy) + X2 cos I ^^i^ ) sin {x2y) 



2k 



. , ^ . , X2W\ . XiW , . , , 
2 sm (xiy) sm ( ] - sm ( ] sm {x2y) 



Akx 



Akx 



(49) 



(50) 



Uy {x,y) 



-A, 



2X1 



cos ixiy) sin 



X2W 



X2 sm 



XiW 



COS (X22/) 



Akx 



(51) 



where xi = ('^^/Cf ~ k^)^^"^ and X2 — (^^/cf — /c^)^^^, with Cf the transverse sound velocity and q the longitudinal 
velocity in a large thin plate 



E 



E 



Ct 



2/9(1 + a)' 



c/ 



p(l-a2) 



(52) 



and uj and k related by the dispersion curve which must be found numerically. In the thin plate limit it is sufficient 
to take for the in-plane modes 



Ux{x,y,z) ~ Ux{x,y), 
Uy{x,y,z) ~ Uy{x,y), 
Uz{x,y,z) ~ 0. 

Similarly, the vertical displacement field for the even flexural modes is 



Uz ix,y) = A2 



cosh 



X-W 



^ cosh (x+y) - _ cosh (^^^^^ cosh (x-y) 



Akx 



(53) 
(54) 
(55) 



(56) 
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and for the odd flexural modes 



Uz (x, y) = A2 



sinh 



X-W 



sinh (x+y) 



Pa 



k'^a — X 



■ sinh 



X+W 



sinh ix-y) 



Akx 



(57) 



where x+ = (k^ + \/pd/Du?j and x~ = (k^ — \/pd/Du?j , with D = EdP / 12 (l — tr-^) the flexural rigidify, and 
again w and k are related by the appropriate dispersion curve. In the classical thin plate theory, the displacement 
fields are given in terms of Uz by the expressions 



Ux{x,y,z) ~ ~zdxUz{x,y), 
Uy{x,y,z) ~ ~zdyUz{x,y), 
Uz{x,y,z) ~ Uz{x,y). 



(58) 
(59) 
(60) 



This approximation is adequate for evaluating the surface stress integrals in Eq. (|4g) but turns out not to be suffi ciently 
accurate to evaluate the energy flux expression for the group velocity Eq. (|2^). We discuss this case in section III C 2 
below. 



A. Ideal thermal conductance 



Since our quantitative calculation of the scattering coefficient relies on the analytic expressions for the elastic modes 
available only in the thin plate limit, it is essential to estimate the temperature range where the thin plate limit is 
applicable for a given experimental structure. On the other hand, as the wavelength becomes much smaller than the 
dimensions of the structure, we should to be able to treat the waves in terms of separate longitudinal and transverse 
waves in the bulk of the material, without worrying too much about the complicated standing wave transverse mode 
structure important for the long wavelength modes. In this regime, which we refer to as the bulk mode limit, the 
counting of the modes is insensitive to the details of the boundary conditions, and is the same as for a scalar wave 
approximation. The ideal thermal conductance depends only on cut-off frequency of the modes (see Eq. (|^)), and 
we can assess the applicability of these simple limiting approximations by comparing the mode cutoff frequencies 
with results-|from a numerical calculation of the full elasticity theory. For the full elastic theory, we use the "xyz 
algorithm"E3. 
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For the bulk mode calculation, there are three polarizations (one longitudinal and two transverse) with propagation 
velocities c^i and Ct, respectively with Ct as before and 



E(l-a) 

C3l ^\ i,_r.t, ' y 61 
y p(l + cr) (1 - 2cr) 

The precise details of the boundary conditions are unimportant in the mode counting for large mode numbers. If we 
assume standing waves in the transverse direction corresponding to zero normal derivative boundary conditions on 
the wave functions, the cutoff frequencies are 



+ ' (*^° degenerate) (62) 



for the transverse waves, and 



ITITT \ / TITT \ 

uJi,mn = C3i\l (^"^j + [-^ ) , (nou degenerate) (63) 



for the longitudinal waves, with m,n = 0,l,2.... For large m, n we can use the continuous form for the frequency 
ujjy oi the A^th mode 

dW ^ f 2 l\ 



Figure || shows the cutoff frequencies as a function of mode number for the thickness to width ratio d/W = 0.38. 
The thin plate theory gives a good approximation at lower frequencies. The accuracy of thin plate theory becomes 
better as d/W gets smaller. For example, in the case of d/W = 0.1 (not shown), the error in the cutoff frequencies of 
the first 13 modes is less than 3%, whilst the error is as large as 5% for the first 7 modes for the case d/W = 0.38 shown 
in the figure. In terms of the ideal (no-scattering) thermal conductance (Eq. (|l|) with the transmission coefficient 
set to unity), we find that for d/W = 0.38 the error in the thermal conductance is less than 4% up to T ~ 0.4 K. 
Thus the thin plate limit is adequate to examine the scattering effects in this temperature range. At large frequencies 
LuW/ct > 30 the elasticity theory results approach closely the continuum bulk mode calculations, Eq. (|64|). The thin 
plate approximation clearly fails in this limit, since it predicts N on uj corresponding to a 2d structure. 

B. Attenuation coefficient in the thin plate limit 

The thin plate approximation is implemented by noticing that the stress free boundary conditions imply that the 
stress components Tiz are zero on the top and bottom surfaces. For small thickness this implies that,the components 
Tiz for any i are small everywhere. In most situations these components can be approximated as zerdij. This simplifies 
many of the terms appearing in Eq. (^l]). Also, at low temperatures, only modes with no strong dependence on the 
z coordinate will be excited, so that the mode sum extends over modes with increasing numbers of nodes in the y 
direction only. 

In this section we calculate the scattering of the elastic waves by surface roughness for a thin plate. We assume 
that the roughness is confined to the sides, since in the experiments theses are prepared lithographically, whereas the 
top and bottom surfaces are produced by the epitaxial growth process. 

For simplicity, we assume the roughness function /i has no z dependence. This is probably a reasonable description 
of the roughness produced by a typical lithographic process of anisotropic chemical etcho. Then the Fourier trans- 
formed roughness function fi {km — kn) can be pulled outside of the z integral in Eq. ( ^ ) and the statistical average 
over the roughness can be performed to give 

h{k)\^^~g{k)L, (65) 
where g (k) is the Fourier transform of the roughness correlation function 

~g{k)^ f dxe-'^-(/i(:E)/i(0)). 
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Equation (Efl) leads to the back-scattering rate from mode m to mode n 



In 



{m) in) 



dzr(™'")(z) 



(66) 



where Eq. ( ^6|) is multiphcd by a factor of 2 to account for the two surfaces B.i y — zkW/2. 

With the closed forms of the displacement fields at hand, we can obtain the analytical expression for the attenuation 
coefficient. We first evaluate f from Eq. (E3). Since t;^") ~ 0, the expression for T reduces to 



V=W/2 



In addition, putting Ti™-* in Eq. ( ^ ) at the stress free boundary to zero gives 
so that Ti™-* from Eq. (^5|) simplifies to 



Now Eq. (|66D can be written as 



g (km kn) 



(67) 



(68) 



(69) 



(70) 



where the index i is summed over x, y, z. The scattering in the thin plate limit is seen to have two components: the 
kinetic term, the first term in the [] in Eq. ([zoj), which involves all components of the displacement; and the stress 
term, the second term, which just depends on the longitudinal displacement. 

To see how the scattering rate scales with the parameters it is useful to rewrite Eq. ( [70|) as 



g{km-K)L W^uj^ 

-fn,mL = TTTT-. X 



2M^4 



/- 



d/2 
■d/2 



dz 



/.d/2 , cW/2 ,(m)*y/2 / .d/2 , fW/2 dy An) An)*V''^' 

\J-d/2"-'^ J-W/2 W^i fi ) \J-d/2^^ J-W/2 w'Pt fi ) 



(71) 



The first factor is a dimensionless measure of the strength of the roughness; the second factor is a dimensionless ratio 
that depends, through the dispersion relation, only on the geometric ratio d/W and the Poisson ratio cr; and the final 
factor involves integrals over the displacement fields, where we have introduced the explicit normalization factors in 
the denominator so that we may evaluate the ratio using convenient unnormalized expressions for the displacements. 



C. Evaluating the group velocity 



As we have seen in Eq. (^0|), we can avoid evaluating the group velocity appearing in Eq. ([70|) via numerical 
differentiating the dispersion curve by instead relating the group velocity to the energy flux in the mode, which in 
turn can be written as an explicit integral Eq. (|2^). Thus we need to evaluate the expression (we suppress the mode 
index in this section ) 



P 



1 



■Re 



dydz 



(72) 



involving the displacement fields and their derivatives. 

In the thin plate limit the z components of the stress are small. If we approximate Tzz = then expressions Eqs. 
(^ and (|^) can be used to evaluate the z-component of the strain 



dzU^ = - 



(1-a) 



{dxUx + dyUy) . 



(73) 
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This can be then used to simplify the in-plane components of the stress 

;{dxUx + adyUy), (74) 



E 


(1 






E 


(1 






E 



Tyy = _ ^2^ i'7^xU^ + dyUy), (75) 
E 



These expression are used to evaluate the first two terms in the integrand in Eq. ([72|). The evaluation of the last 
term in the integrand turns out to depend on whether we are looking at the in-plane or flexural modes, and we now 
consider each case in turn. 

1. In-plane modes 

For the in-plane modes in the thin plate limit it is sufficiently accurate to approximate T^x — 0, and we can evaluate 
the remaining terms in P with the approximations — u^^ Uy Uy independent of z. This yields 

P = Re 1^ J^Y~~5y j '^y [^{dxUx + crdyUy)ul + (1 - cr) {d^Uy + dyUx)ul] | . (77) 

2. Flexural modes 

For the flexural mode the approximations Tzx — and Uz{x, y, z) ~ Uzix, y) independent of z lead to the expressions 
for the horizontal displacements 

Ux{x, y, z) ~ ~zuz{x, y), (78) 
Uy{x, y, z) ~ -zuzix, y). (79) 

Using these expressions with Eqs. (|7^[75|) shows that the first two terms in Eq. ( |7^ ) are of order d^, i.e. third order in 
the expansion parameter of thin plate theory d/W . It turns out that to this order, we cannot neglect the last term in 
Tzx-, even though all z-components in the stress tensor are nominally "small". Indeed comparing the group velocity 
evaluated from Eq. (^) neglecting the term in Tzx with those given by numerically differentiating the dispersion curve 
shows a clear discrepancy. This same problem comes up in deriving the wave equation for the flexural waves 

pdui^Uz = D\/iuz. (80) 

The term on the left hand side is the mass per unit area times the vertical acceleration, which is given by the integral 
over the depth of dxTzx +dyTzy. Clearly the components of T^^ cannot be neglected completely. Their "smallness" is 
what leads to the unusual fourth order derivative appearing in this wave equation, with a coefficient again proportional 
to d^. 

We have used two methods to arrive at the correct calculation of the energy flux integral for the flexural waves, which 
is then used to calculate the group velocity for the these waves. The first is to use an improved approximation to the 
expressions for the in-plane displacements Eqs. ( |7^ , |75| ) and a nonzero Tzx following the approach of Timoshenkocil. 
The second evaluates the energy flux in terms of the vertical displacement and an effective vertical force--and in 
addition the rotational displacement 9 and corresponding torque M, as is used in the macroscopic derivationt^ of the 
wave equation (^) . Either of these methods leads to the expression for the energy flux 



P ~ -LuD Rc 
2 



^ivDRc I^J dy [2k^UzU*z + fc(l - a){dyUz)dyU*z - fc(l + (T){d^Uz)u*z] 

+Dk [(1 - cT){dyUz)ul]y^^ - Dk [(1 - a){dyUz)u*z]y^^^'} (81) 



The derivations are displayed in Appendix The comparison of the group velocity derived from Eq. ( |8l| ) and from 
numerically differentiating the dispersion curve now shows agreement to high accuracy. 
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FIG. 3: Dispersion relation for in-plane modes (solid) and flexural modes (dashed) for a geometry ratio d/W = 0.375 and 
Poisson ration 0.24. The wave numbers are scaled with the width W, and the frequencies by W/ct with ct = \/ p^j p. 



IV. SCATTERING ANALYSIS 



The thermal attenuation is calculated from Eq. (|7y) for normalized mode displacement fields or ( |7l| ) in general. 
The group velocity for each mode can be accurately evaluated numerically from the equality Vg = 2P/ pw^ , with the 
energy flux P given by Eq. (|7^) for the in plane modes and Eq. (|8l|) for the flexural modes (both expressions are for 
normalized displacement fields) . These are all explicit results in terms of the mode displacements, which are given by 
Eqs. (^8[]5l| ) for the in-plane modes, and Eqs. ( p6| , ^7|) for the flexural modes. 

Before analyzing the scattering behavior, we first need to have a good understanding of the dispersion relation of 
the modes, since the scattering rates are strongly dependent on this. 



A. Dispersion relation and group velocity 

The dispersion relations for a representative case are shown in Fig. (^. For this example we have used a Poisson 
ratio of 0.24, and a depth to width ratio of d/W = 0.375, values corresponding to the experimental work of Schwab et 
al.H. As we have discussed, the modes fall into four classes, depending on their parity signatures. We label the lowest 
mode from each class, the one with zero frequency as the wave number goes to zero, as mode 0, and the modes with 
successively higher cutoff frequencies in each class as mode 1, mode 2, etc., in that class. 

Notice that one of the curves in the figure, the one for the in-plane mode with cutoff frequency uiW/ct — 5, 
shows anomalous dispersion with the frequency decreasing as the wave number increases up to about 3M^~^. (This 
is actually an even mode, and some higher even and odd modes also show anomalous dispersion.) The dispersion 
curves for all modes n > have zero slope, and so zero group velocity, at onset. As we will see later, this results in 
a diverging scattering rate at each mode onset. For the n — Q modes, as w — s- two of the modes (the compression 
and torsion modes) have linear dispersion, whilst the other two lowest modes (in-plane and flexural bending modes) 
exhibit quadratic dispersion. Figure (|^) shows the group velocities Vg for the four lowest in-plane modes. The group 
velocity of the bending mode approaches zero as w — > whilst that of the compression mode becomes constant. The 
group velocity of the compression mode suddenly drops to ^ O.Scj around ojW/ Ct 4.6, then gradually recovers and 
approaches 0.9ct. These features of the dispersion curve will be reflected in the behavior of the scattering of the waves. 



B. Scattering behavior 

We first consider the scattering and reduction of the thermal transport by white noise roughness g{k) = 5(0). This 
allows us to focus on the role of geometry and the unusual mode structure of the elastic waves in the physics. 
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FIG. 4; Group velocity for in-plane modes for the same parameters as Fig. (^): dash dotted - in-plane bending mode; solid 
compression mode. The wave numbers are scaled with the width W, and the group velocities by Ct with ct = \J i^j p- 



Extension k 11 0{ky) 0(kz) 

In-planc bend {w / ^/V2)k'^ {w/V3)k -iky 1 0{kz) 

Torsion ^2/(1 + a)(d/w)k ^2/(1 -|- (7)(d/w) 0{kyz) -z y 

Flex-bend (d/VT2)fe^ (d/V3)k -ikz O(k^yz) 1 



TABLE L Dispersion relation, group velocity, and (unnormalized) transverse mode structure for the four modes with zero 
frequency at zero wave vector. 



In the low frequency limit the dispersion curve and the spatial dependence of the modes take on the simple analytic 
forms shown in Table |, allowing us to make analytic predictions for the scattering at low frequencies, and then the 
thermal conductance at low temperatures. Since only small wave vector scattering is involved in these calculations, 
the results are true for a general roughness correlation function, providing g{0) is nonzero. The mode structure in 
Table ^ may be calculated from Eqs. (^^57) taking A; or from arguments of macroscopic elasticity theory. 



The contributions to the thermal attenuatiojt-cocfficicnt in the low frequency limit {coW/ct <C 1) from the various 
scattering processes are shown in the Table pH3. The expressions take on their simplest form if we introd uce t he 
frequency scaled with the velocity of the long wavelength compression mode uj = luce/W with ce = \/ E/p — 
•\/2(l -|- a)ct. The power laws can largely be understood from the prefactor in Eq. (|7l]), jn,m oc 

The 

group velocity Vg becomes a constant at small frequencies for the compression and torsion modes. Thus the torsion- 
torsion and compression-compression scattering shows the dependence corresponding to Rayleigh scattering in 
one dimension, and as was found for scalar waves with linear dispersion. On the other hand for the bending modes 
Vg cx w^/^. This has the important consequence that the in- plane bend-bend scattering increases more rapidly at 



in-plane 



bb 



23/2'- 



flexural 



0(1+J) t W: 



bb 



O {1 



P747 



tb,bt 



TT7T 



TABLE II: Scattering coefficients for the zero onset frequency modes at low frequencies: c denotes compression, b denotes bend, 
t denotes torsion, bb denotes bend to bend scattering etc. Values are quoted for 7mW^*/ff(0) as a function of scaled frequency 
uj = ujce/W. For the flexural bend to bend scattering (bb) the terms in the braces in Eq. (hu) cancel to leading order resulting 
in very small 0(a;^) scattering. There is no scattering between in-plane and flexural modes for the z-independent roughness 
assumed. 
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FIG. 5: Attenuation coefficient jmW* /g{0) for scattering from the two lowest m = inplane modes to any other mode as a 
function of scaled frequency coW/ct: solid line - inplane bend mode; dashed line - compression mode. The insert shows an 
enlargement of the low frequency region, and compares with the analytic low frequency expressions from Table [l[ dotted - 
analytic inplane bend mode; dash-dotted - analytic compression mode; other lines as in the main figure. 



low frequencies proportional to w, and the torsion-bend and compression-bend scattering have an w^^^ frequency 
dependence. For the flexural bend-bend scattering the two terms in the braces in Eq. d7G) cancel to leading order 
resulting in smaller scattering 0{uj^) than given by the prefactor alone. Note that the expressions for the flexural 
modes involve additional factors of W/d, so that these modes will be scattered more strongly at a given uj in the 
thin plate limit. This is because these modes are softer, so that the scattering wave vectors are larger for the same 
frequency. 

Numerical results for the attenuation coefficient 7™ of the four lowest modes are shown in Fig. (^) for the in-plane 
and Fig. for the flexural modes. The plot for the in-plane modes in particular shows interesting structure deriving 
from the complicated dispersion curves of Fig. ^. Much of this structure can be understood from the product of 
group velocities in the denominator of Eq. (|7^). In particular there is a square root divergence in 7™ at the onset 
frequency of each mode where the group velocity is zero. In addition, the large scattering around LoW/ct — 5 derives 
from the region of anomalous dispersion, since the group velocity is small in this frequency range. The insert to Fig. 

shows an expanded view of the low frequency behavior, using the results from Tabic |l| together with the next 
order correction for the compression-bend scattering. The agreement for the compression mode is very good even up 
to iuW/ct ~ 3, whereas for the bend mode the correspondence is only good for ujW/ct ^ 0.5. The scattering for the 
flexural modes shows generally similar results. Fig. (g) although the behavior is simpler corresponding to the rather 
featureless dispersion curves. At low frequencies, (insert to Fig. (|^)), the scattering of the flexural-bend mode is small, 
since the intramode scattering is reduced by the cancellation discussed above. 

Figure (j^) shows the total scattering 7^ for the in-plane modes on a log-log plot, again with white noise 
roughness. At very low frequencies the scattering varies proportional to co corresponding to the dominant intramode 
scattering of the compression mode at low frequencies (Table ||). For frequencies up to ujW/ct — 3.5, the flrst 
nonzero onset frequency of an in-plane mode, the analytic low- frequency expression given by summing the in-scattering 
expressions from Table || (cc, cb, be, and bb), shown as the dotted line in Fig. |^, gives a good approximation to the 
full results. At higher frequencies the total scattering increases rapidly, following a general trend proportional to 
(dashed line) together with divergent scattering at each mode onset frequency. The uj'^ power law can be understood 
as the combination of the explicit dependence of Eq. ( [70[ ) , together with two powers of uj coming from the number 
of modes available for scattering from and to. 
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FIG. 6: Attenuation coefficient JmW' /g{0) for scattering from tiie two lowest m = flex modes to any other mode as a function 
of scaled frequency lu^12{1 — a'^){W / d)W / ce- solid line - the flex-bend mode; dashed line - torsion mode. The insert shows 
an enlargement of the low frequency region, and compares with the analytic low frequency expressions from Table dotted 
line - analytic approximation for the flex-bend mode; dash-dotted: analytic expression for the torsion mode; other lines as in 
the main figure. 



C. Change in the thermal conductance 



In the weak scattering limit the change in thermal conductance at low temperatures can be derived directly from 
the expressions for the scattering at low frequencies. If we write the thermal attenuation coefficient of mode m as 
7„ii — A{lo/woY, where p is the power law obtained in the low frequency limit and some characteristic frequency, 
then the corresponding contribution of the suppression of the thermal conductance from this mode is 

SKjKu = AIp{T/T^Y (82) 

with Tq — htoo/kB the corresponding characteristic temperature and K^^ = Tr^kgT/3h the universal thermal conduc- 
tance. The constant Ip can be obtained evaluating the integral 

Ip = — dy— —. (83) 

Jo ye-y - ly 

Thus the power law for the temperature dependence of the depression of the thermal conductivity is the same as the 
one for the low frequency behavior of the scattering coefficient. 

Figures (||) and ^) show the thermal conductance depression scaled with the universal value as a function of 
the appropriate scaled temperature for the lowest in-plane and flexural modes, showing the deviation from the low 
temperature power laws as the temperature is raised. For the in plane modes we use the characteristic temperature 
Te = hcE/ksW and for the flexural modes Tp = hcEd/ksW^ . The individual plots are then independent of the 
geometry. To combine the contributions from the in-plane and flexural modes the ratio d/W \s needed to relate the 
two temperature scale factors. In the thin plate limit Tp = {d/W)TE <C Te- 



V. COMPARISON WITH EXPERIMENT 



A. Experimental geometry 

Based on the SEM micrograph of the experimental structureS, we set the dimensions of the structure in the 
following way. In the experimental structure of Schwab et al. the thermal pathway was constructed with the shape 
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FIG. 7: Total scattering 7mH^*/5(0) for the in-plane modes on a log-log plot. The dotted line shows the low frequency 
analytic expression from Table ^ and the dashed line shows a power law 4. (Note that the heights of the peaks in the plot are 
not significant, depending on how close the individual points, separated by 0.01 in uW/ct, used in constructing the plot are to 
the mode onset frequencies where the scattering diverges.) 
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Scaled temperature T/T^ 
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FIG. 8: Reduction in the thermal conductance scaled with the universal conductance Ku for the lowest in-plane modes as a 
function of scaled temperature T/Te with Te = hcE/ksW: solid line - low temperature analytical expressions from Table ^ 
points - full expression evaluated numerically. The quantity plotted is {SKc + 5Kii,)/2Ku with SKc, SKu, the depression of the 
contributions to the conductance by the scattering for the compression and in-plane bending modes. 



function W{x) = W cosh(Ax) so that the beam width becomes large and joins smoothly to the thermal reservoirs 
at the ends, reducing the scattering due to the geometric imperfection at these junctions. Unfortunately this makes 
the calculation of the behavior of the elastic waves in the beams much harder. However both with and without the 
scattering off surface roughness, we expect the narrow portion of the beam to dominate the behavior. Thus we simplify 
the structure and model it as an elastic beam with rectangular cross section of width W , depth d, and effective length 
L. We estimate the width as the narrowest width of the structure, W ~ 160nm, and L = 1 ^m as the length over 
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Ssaled temperature (W/d)T/Tg 

FIG. 9: Similar to Fig. SK/2Ku for the lowest flexural modes (torsion and flexural-bending) as a function of the scaled 
temperature T/Tf with Tp = HcEd/ksW^. 



which the width is approximately constant. The thickness of the material was d — 60 nm. The accuracy of the length 
estimation is not very critical, since the only length dependence in the scattering rate 7 appears in the combination 
S^L where 5 is the rms roughness which is a parameter of the model, so that any error in the assignment of L will 
just change the value assigned to d. The width W on the other hand plays a crucial role, for example determining 
the frequency cutoffs of the various modes, and so the temperature dependence of the thermal conductivity. 



B. Roughness correlation function 



Since the nature of the surface roughness on the experimental structure is not known, to fit the experimental data 
we need a sensible parameterization of the roughness. As a starting point we choose a Gaussian correlation function 
for the roughness, leading to the spectral density 



g (fc) = y/naS'^ exp 



(84) 



This parameterization of the roughness contains two parameters: S the rms roughness and a the correlation length. 

To analyze the data, we first quantify the amount of scattering by subtracting the data of Schwab et al. from the 
ideal thermal conductance obtained numerically using the "xyz" algorithmEJ. Then we attempt to fit the data by 
adjusting the two parameters a and S^L. 

The inadequacy of Eq. (|^) in fitting the experimental data is shown by the low temperature fits in Fig. (10). At 
these low temperatures only small wave number modes are excited, so that the exponential term in Eq. (p3) can be 
approximated as unity and g (k) ~ g (0) — ^piiaS^ . Thus the roughness parameters only appear in the combination a(5^, 
and this quantity can be varied to attempt to fit the low temperature region. As seen from the figure, increasing aS^ 
causes scattering that is systematically larger than the experimental data at the low temperatures, while decreasing 
aS^ does not provide enough scattering in the range 0.1 < T < 0.2 K. 

Although there is considerable scatter in the data over the range of the fit, the systematic differences between the 
predictions and the data lead us to propose a modified form of the roughness correlation that reduces the scattering 
at small wave numbers 



g (k) = ^/naS'^ exp 



(k - koY 



(85) 



A nonzero value of the parameter fcp leads to a roughness correlation function that is maximum at a length scale of 
order k^^ , and serves to reduce the scattering at long wavelengths. As mentioned in the introduction section, the same 
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FIG. 10: Attempts to fit the low temperature data T < 0.2 K using various values of aS^: solid line - y'TiaS^ — 0.1; dotted line 
- V7ro(5^ = 0.05; dashed line - y/^aS'^ = 0.02; open circles - from the experimental data of Schwab et al. 
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FIG. 11: Thermal conductance per mode scaled with universal value Ku'- solid line - fit using roughness parameters a/W = 5.5, 
d/W = 0.2, and koW = 4.9; circles - data of Schwab et al. The dotted line shows the ideal value with no scattering. 



discrepancy (i.e., the overestimation of the scattering aJ^long wavelengths in the theory compared with experiment) 
was found using the scalar model of the elastic wavesllj. The full elasticity theory considered here actually makes 
the discrepancy worse, since the scattering at small frequencies now is predicted to increase more rapidly at small 
frequencies than the found in the scalar theory, varying as ujP with p < 2 for most of the scattering processes, see 
Table || 

To fit the data of Schwab et al, we need to determine three parameters: fcg, a, and S. As well as looking at the fit 
by eye, we evaluate the quality of the fit by calculating the mean square deviation of the data from the theory curve 
over the temperature range up to OAK . At higher temperatures many modes becoming excited, and the scattering 
of individual modes becomes strong, so that our theory is less reliable. Since the onset frequency of the scattering at 
low frequencies and the initial decrease in thermal conductance with increasing temperature near the onset is mainly 
determined by fco, this parameter is the easiest to determine. We find the value k^W = 4.9 , rather insensitive to the 
values of a and S. We have also investigated the fit with a delta function noise correlation function (i.e. a — s- in 
Eq. (|8^)). With this roughness, to get a maximum in the decrease in the thermal conductance at about T ^ 0.2K 
requires kgW values between 4.5 and 5, which agrees with the previous estimate. 

The remaining fit parameters a and S are not well determined. Visual inspection of the fit and a plot of the error as 
a function of {a,S) shows a reasonable fit over the range a/W = 5,S/W — 0.17 to a/W = 6.7, 6/W = 0.32. A choice 
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FIG. 12: Same as in 



Fig. (hll) but showing the decrease of K/Ku from the ideal value. 



of the roughness amplitude towards the lower end of this range seems most likely physically (e.g., S/W = 0.31 ~ 50 
nm is unlikely for roughness generated by chemical etch). We therefore use S/W = 0.2 and a/W — 5.5 to generate 
the curves. 

A difficulty of fitting the data is the lack of data points at very low temperatures: it is in this range where 
only a few modes are involved that we have a very good understanding of the scattering. At higher temperatures 
many more modes become involved, and the scattering of individual modes becomes strong, so that the second order 
approximation used in calculating the scattering will not be good. A full test of the theory explaining the reduction 
in the thermal conductance in terms of the scattering off surface roughness requires more data below a temperature 
of about 0.Q8K for the type of geometry used by Schwab et al., or systems with smaller geometries where the effects 
can be measured at higher temperatures. 



It is interesting to investigate the contribution to the total thermal conductance of the individual modes with the 
roughness parameters used to fit the experimental data. This is shown in Fig. ([l3|). The flex-bending mode shows 
a much smaller contribution to the reduction in K at low temperatures for the reason we have already discussed. 
The modes with nonzero onset frequencies start to contribute significantly above about T ~ 0.2K, and this is the 
predominant cause for the increase in thermal conductivity above this temperature, since the recovery of the thermal 
conductance for the lowest mode occurs very slowly. 



We have investigated the effect of surface roughness on the scattering of elastic waves in a rectangular beam or 
waveguide, and the resulting depression of the thermal conductance in the low temperature quantized limit, using full 
elasticity theory. Our formulation is quite general, but to obtain concrete results we have specialized to the thin-plate 
limit, which should be a reasonable approximation for many mesoscopic experiments where the depth of the structures 
is fixed by the epitaxial growth, whilst the width is determined lithographically. The thin plate limit preserves the 
peculiar features of the elastic waves in the full elastic theory, namely a quadratic dispersion at long wavelengths 
for two of the low frequency modes, and regions of negative dispersion in the spectra. A robust result is that the 
low frequency asymptotic dependence of the scattering by unstructured roughness of the modes that propagate at 
low frequencies (the ones that are important in the low temperature universal thermal conductance) depends on the 
structure of the modes and the dispersion relation, and is not the simple dependence of Rayleigh scattering as found 
in the scalar approximation to the modes. We find different power laws for the various mode scattering processes that 



C. 



Individual mode contribution to the thermal conductance 



VI. CONCLUSION 



22 




0.2 - 



0.0_ 



0.0 



0.1 



0.2 



0.3 



0.4 



0.5 



T(K) 



FIG. 13: Individual mode contribution to the thermal conductance. The lowest two flex modes and lowest three in plane modes 
are shown. The contributions to K/Ku from the four modes with zero onset frequency tend to unity at low temperatures. The 
higher modes only contribute at higher temperature. The modes are: dash-dotted - in-plane bending; dashed - compression; 
dotted - torsion; dashed-dotted-dot - out-of plane bending. The solid line shows the sum of all the mode contributions, reduced 
by AKri. Values of the roughness parameters used were a/W = 5.5, 5/W = 0.2, koW — 4.9, and d/W — 0.375. 



can be understood largely from the dispersion relations at: lo for intramode scattering for the in-plane bend mode 
(the flex-bend intramode scattering is anomalous because of a cancellation between leading order terms, and varies 
as w'^); u;^/^ for scattering between the bend modes and the modes with linear dispersion (torsion and compression 
modes); and the usual for the intramode scattering of the modes with linear dispersion. The cm'rent experimental 
data on the suppression of the low temperature thermal conductance below the universal value does not extend to low 
enough temperatures to provide a good test of these predictions. To investigate this prediction further, it would be 
interesting to extend the experiments to lower temperatures, or to smaller devices such as carbon nanotubes, where 
the characteristic temperature scales (when a typical thermally excited phonon has a wavelength comparable with 
the device dimensions) are higher. 

We have used our results to understand the data the data of Schwab et al., who observed a depression of the 
thermal conductance below the universal value in the temperature range of O.IK to OAK. Although the scatter in the 
data is considerable at these low temperatures, the observations seem to show a delay in the onset of the depression 
scattering as the temperature is raised, beyond what can be fitted with our predictions for unstructured surface 
roughness. We tentatively resolve this delay by supposing that the surface roughness has a maximum amplitude at 
some nonzero length scale, which we parameterize by a shifted Gaussian correlation function. Due to the lack of data 
at low temperatures, a precise determination of the roughness parameters is not possible. However, we do obtain a fit 
to the data with parameters that do not look unreasonable when compared with electron micrographs of the actual 
devices. 

Our results are based on second order perturbation theory, and the thermal conductance is evaluated assuming the 
scattering over the length of the device is small. This is a good approximation at low temperatures, but the scattering 
becomes strong at higher temperatures, particularly for the new modes excited as the temperature is raised, which have 
a diverging scattering at onset due to the flat dispersion relation here. At higher rtemperatures multiple scattering 
and perhaps phonon localization will therefore become important. Kambili et al.ll3 and Sanchez-Gil et al.tll have 
numerically investigated the these effects in the simplified scalar wave approximation. It would be interesting in the 
future to extend their work to the full elasticity model. 
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APPENDIX A: INCIDENT AND SCATTERED FIELDS 



Using Green's theorem we have expressed the displacement field at frequency uj in terms of the surface integral 

M,(x) = / [n'r„- (x') (x', x) - n'.u, (x') Ty, (x', x)] dS' (Al) 

JS' 



Eq. ( Al) involves the integration over a closed surface S", which we have chosen to be the smooth boundaries together 
with the cross sections at x' — s- ±cx). In this appendix we show that the integration over the sections at ±00 simply 
yields the incident field w™, and this allows us to deduce the expression for the scattered field as an integration over 
the side, surfaces. To deduce this result, we first need to derive what arc known as reciprocity relations for the elastic 
modealj. 

Let u^'') and u^^^ be the displacement fields for modes r and s in the ideal beam, and T^''^ T^'*) the corresponding 
stress tensor fields. The modes satisfy the wave equation at frequency uj, so that 



9 (r) 
2 (s) 



(A2) 



Multiply the first equation by u[^^* and the complex conjugate of the second by u'f\ subtract the two equations, 
integrate over a volume of the beam between x — xi and x = X2, and finally use the divergence theorem to find 

fijdS = 0, (A3) 



(r)r7i(s)* 



where the integral is over the surface bounding the volume, consisting of the sides of the beam between xi and X2^ and 
the sections at Xi and X2- The integrations over the sides of the beam are zero by the stress free boundary conditions. 
For the integration over the sections introduce the explicit x-dependence u^*"^ — (j){y, z)e^''''^ and T^''^ = T^'^\y, z)e^^''^ 
with kr the wave number of mode r at frequency u) etc. Then Eq. ([A3|) reduces to 



dydz = 



(A4) 



and the integral is independent of x. Unless the prefactor is zero, this shows us that the integral over the section must 
be zero, and so 



dydz = 0, kr ^ kg 



(A5) 



This is one version of the reciprocity relations. 

For our purposes it is more convenient to express the condition for the reciprocity integral to be zero in terms of 
the group velocity rather than the wave number. To do so, we need to consider the dispersion curves. The condition 
for the reciprocity integral to be nonzero, k^ = kg for modes r, s at the same frequency w, actually implies r and s 

(r) (s) 

are the same mode, so that in fact Vg — Vg . The only other possibility is that r and s are modes with dispersion 
curves that cross at frequency uj, k = kr = kg. However only modes of different y, z parity signatures can cross, and 
then the integration over the section for these different modes in Eq. ( |A5D is again zero. Thus we can rewrite the 
reciprocity relation as 



{s)*rr,(r) 



dydz^O, v'^pT^v'fK (A6) 
If r and s are the same mode, the integral is related to the energy fiux and hence to the group velocity (see Eq. dlf 



j j dydz 



{r)*rr(r) 



2ipujvpK 



(A7) 



We now use Eqs. (A6,A7) to evaluate the contributions to Eq. (Al) from the integrations over the sections at 
x' — > ±00. 

Let us first consider x' 00. According to Eq. (|l8|) the x' dependence of the Green's function pair G, F consist of 
modes Us(a;')* with Wg'*-' < since here x' > x for any finite x. On the other hand the field pair u, T are made up of the 
incident wave, and waves scattered from the roughness at finite x, and so consist of modes Ur{x') with Ug'^^ > 0. The 
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integral in Eq. ( Al) over the section at — > oo is therefore the sum of terms involving J J 



with Vg^^ and vf of opposite sign. All these terms are zero by Eq. (A6), and so there is no contribution from the 
section at a;' ^ oo. 

Similar arguments apply to the section at x' — oo. The Green function is made up of modes with Vg > 0. The 
scattered component of the field u consists of modes with Vg < 0, and there is no contribution to the integral over 

the section from these modes. On the other hand the incident wave u'" is mode Um with ui™'' > 0, and there is the 



„(«) 



Ur 1> 



dydz 



(n) 



single term with v, 
gives u^q^\x). Writing u 



surviving in the sum over modes in the Green function. Using Eq. (A7) the integral just 
u'" + u"'^ then leads to Eq. (|l|) in the text. 



APPENDIX B: ENERGY FLUX FOR FLEXURAL MODES 

The classical thin plate approximation of setting Tzi = is not sufficient to calculate the energy fiux of the flexural 
modes using the integral Eq. ([72|). In this appendix we evaluate the correct-expression for the energy flux by two 
different methods, first using the extended thin-plate, theory of TimoshenkcO (see also GrafS), and then using a 
method in terms of the energy of plate deformationsE3 that avoids these difficulties. 

In the extended thin plate approximation of Timoshenko the z-dependence of the in-plane displacements is still 
approximated as linear 

Ux{x, y, z) ~ z^x{x, y), (Bl) 
Uy{x, y, z) ~ zi/jyix, y). (B2) 

However, the x, y dependence is no longer assumed to be given by the gradient of the mean vertical displacement 
Uz(x,y), but by the more general expression 

^ = -Viu^+V±5 + V± X (Cz) (B3) 

introducing the scalar potential S{x^ y) and vector potential C^{x^ y) defining the corrections to the in-plane strain and 
rotation. Here Vj_ — {dx, dy) is the horizontal gradient. In addition the vertically averaged stress Tzx is taken to be 

Tzx ^ n'—^—idxUz + dzu,) (B4) 

(with a similar expression for Tzy given by replacing the subscript x with y everywhere). Here the "shear correction 
factor" K, a number of order unity, is, introduced to take into account deviations of the in-plane displacements from 
the assumed linear dependence on zEd. In the usual thin plate approximation Tzi are set to zero and i}} = — V_l?x;, so 
that {ux,Uy) = — zV_lw. 

With the Timoshenko approximations, the equations of motion for the three components of displacement are now 
investigated. j— ■ 
The equations of motion for the horizontal displacement lead to an equation relating to ujE^ 

y {(1 -cr)VV + (1 + o-)Vj.Vj_ • V'} - ^^^^J^d{iJ + \/ ^Uz) = (B5) 
(remember D ~ Ed^ /\2{\ — cr^), with E Young's modulus, and the shear modulus). The inertial terms d^ip turn 



out to be neghgible in this equation. Using Eq. (|B3|), Eq. ( |B5| ) becomes 



DV^lVKS - w) - K^^d\/j_S + y (1 - ct)V_l X (ViC-z) - kVc^V_l x {(z) = 0. (B6) 



Taking the vertical curl of Eq. (B6) gives 

^{l-a)\/]^n- K^ndn^O (B7) 

with n — z ■ S/± X tjj ~ — V^^C the rotation. For a wave disturbance e**"'^, this gives an exponential dependence on y, 
e^^y with 
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Since A "'^ ~ d <C W, the rotation will be large only over a boundary layer region with width of order d near the edges 
y = zkW/2, where the solution takes the form 



fl{x, y ~ ±W/2) ~ n{±W/2)e"'''e 



ikx -\\vtW/2\ 



(B9) 



The vector potential ^ has a similar solution, so that last two terms in Eq. (B6) cancel. This leaves for the scalar 
potential S 



which immediately gives 



(BIO) 



(Bll) 



(We are only interested in and so do not need to keep track of the arbitrary gradient-free function that could 

be added to this equation.) |— ■ 
The equation of motion for the vertical displacement igH3 



K^^dV\S — —pduj^Uz 



Together Eqs. (B11,B12) give 



pduj'^{uz — 



D 

K^j.id 



Vizi,) = DViu,. 



(B12) 



(B13) 



This is the usual fourth order wave equation, with a small correction term of order {d/W)'^ (the second term in the 
brackets on the left hand side). Note that solutions to this equation vary on the long scale of order fc~^ or W, and 
not the small scale ^ d, so that to a good approximation we have 



pduf^Uz — DV\_Uz 
S ^ 



K^fid 



(B14a) 
(B14b) 



The first equation is now the standard fourth order wave equation. The second equation for S shows it to be small 
compared with by of order (d/W)'^. 

The boundary conditions at the edges are that all stresses are zero, so that in particular a.t y = ±.W/2 



Substituting Eq. (B3) into this gives 



dzTzy — n^pd{dyUz + ipy) — 0. 



dyS - d^C = 0. 



(B15) 



(B16a) 



Equation (B16a) together with Eq. (B14b) tells us the size of the Q correction, which ai y — ±W/2 takes the value 

D 1 



ax,y = ±W/2) 



K^fid ik 



{dyVlu.) 



y=±W/2 ■ 



(B17) 



This expression can be simplified using the boundary condition Tyy — at y — ±W/2, which from Eq. ( [75| ) and Eqs. 
i^M) gives at y = ±W/2 



so that 



ax,y^±W/2) 



K^fid 



The potential C is only large in the boundary layers near the edges where it takes the form 



C(x,y~±W^/2) 



ikD(l-<7) 

K^fid 



{dyUz)\y^^^f/2 ^ 



-\\yTW/2\ 



(B18) 



(B19) 



(B20) 
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Th us fin ally we have expressions for the horizontal displacement field, Eqs. (B2) and (Bl) together with Eqs. ( B14b| ) 
and ( |B20| ) defining S and C,, and Eq. (B4). These can be used to calculate the additional contribution to the energy 

(|72|). (The corrections to Ux and Uy derived here do not change the contributions 
to the order we require, since these terms are already third order in the small 



flux coming from the T^x term in Ec 
from the first two terms in Eq. 
parameter d/W .) 

We therefore need to evaluate 



Tzxuldydz ~ ^d / dy{dxS + dyC,)u 



(B21) 



Both terms in the integral give contributions at the same order. The first term, coming from the correction to the 
in-plane strain Eq. ( B14b| ), is 



K fid dy{dxS)ul 



D J {dxVlu,)uldy. 



(B22) 



The second term in the integrand is only large in the boundary layer region near the edges and from Eq. ( B19| ) 
evaluates to the edge contributions 



K^rf / dy{dyC)ul = -ikD{l - a) [{dyUz)ul]\lJ^{^^ 



(B23) 



Combining these expressions for Eq. ( B21 ) with Eqs. (|7^) together with (|7^) and (^9|) yields the final expression 

dy [2Pu,ul + (1 - a) {dyU,) {dyU,)* - (1 + a) (a^il,) u;] 

+ [{l-a)[dyu.)ul]l=j^^^]. 

which is identical to Eq. 

An alternative approach to calculate the energy flux is to use the expression for the energy of distortions of the 



An alternative approacn to calculate tne energy flux is to use tne express 
plate evaluated using the lowest order expressions Eqs. (|4|j75l) and (^HflP 



F = -D 
2 



(Viu,) +2(l-a) 



dxdy 



dx^ dy^ 



dxdy. 



(B25) 



It turns out that the higher order corrections discussed above are not needed in this expression, and so we can derive 
the energy flux without these difficulties. The functional derivative of F with respect to Uz yields the vertical force 
per unit area in the interior of the plate, which can be used to derive the fourth order wave equation, as well as 
expressions for the energy flux into the plate across the boundaries. The latter expressions give us the result for the 
energy flux along the beam 



P^-Re 
2 



Mx9l + Vul dy + {F^ully^^/^ +FcU*z\y=-w/2) 



where 

V^-Ddx[dluz + i2-a)dluz] 
is the effective vertical force that couples to the vertical displacement Uz, 

Mx = -D {plu, + adluz) 
is the torque that couples to the angular displacement — duz/dx, and 

F,(y = ±W/2) = ± 2D{1 - ^)dlyUz\^^^^^^ 

is a vertical force loc alized at t he edges o f the plate. 
Substituting Eqs. (p27[jB29D into Eq. (p26|) gives 



(B26) 

(B27) 
(B28) 
(B29) 



p = \r. 



(iujD) 



dy {dluz + (rdluz) {-dxUz)* + j dy[dluz + (2 - cr) dxdluz]ul 

-2{\-a){dxdyUz)ul\y^^l^+2{l-a){dxdyUz)u*z\y^_^f^ } (B30) 
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Evaluating dx = ik, and using integration by parts, we again get Eq. (|8l|) 
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Note that the scattering amplitude normalized by the energy flux t„^m = ^j'^s"' | / ji's'"^ jin.m can be seen from Eq. ( ^ ) to 

explicitly satisfy the reciprocity relationr-Sii,™ = For systems in which energy is conserved, as in our case, reciprocity 

is equivalent to time-reversal invariancell3. 

Without this assumption, we would find slightly different z averages of the roughness function /i involved for the scattering 
of the in-plane modes and of the fiexural modes — in fact a direct average for the inplane modes and an average weighted by 

for the flexural modes. In addition there would now be scattering from in-plane to fiexural modes, and vice versa. 
A more accurate expression for the scattering between the in-plane compression and bending modes is giving by keeping the 
next order term in Cu which is 0{lu^^'^) for this scattering process. To include this correction multiply the expression in the 
table by {1 + s/^ /12^''^ f . 



